IUCAA  GW  Workshop 2026
PN  theory  tutorial - 1

Sayantani  Datta

<<xAct`xTensor`

Define  manifold and metric

DefManifold[Σ,3,{i,j,k,l,α, β}];

DefMetric[1,S[-i,-j],DD,{"|","∂"}]

Define constants and tensors

DefConstantSymbol/@{G,c,m,m1,m2,X1,X2,ν,r,rdot, γ,ω};

massRules={m->m1+m2,X1->m1/m,X2->m2/m,ν->X1 X2};

gammaRule=γ->(G m)/(r (c^2));
revgammaRule=(G m)/(r (c^2))->γ;

(*Relative separation vector*)
DefTensor[x[i],Σ]

(*Unit radial vector*)
DefTensor[n[i],Σ]

(*Tangential basis vector*)
DefTensor[λ[i],Σ];

(*Velocity and acceleration vectors*)
DefTensor[v[i],Σ]
DefTensor[a[i],Σ]

COM  definitions/ EOMs

(* Circular-orbit geometry *)
geometryRules={n[i] n[-i]->1,λ[i] λ[-i]->1,n[i] λ[-i]->0};
CircularOrbitRule={rdot->0};

(*Relative position*)
xRule=x[i]:>r n[i];

(* 0PN orbital frequency*)(*Eq.(3.11) *)
omega0PNSqrule=(ω^2)-> (G m/(r^3));
omega0PNpow4rule=(ω^4)-> (((G m/(r^3)))^2);

(*Velocity vector*)
velocityRule=v[i_]:>r ω λ[i];

(*Acceleration vector*)
accnRule=a[i_]:>-(ω^2) x[i];

a0PN[i_] = a[i]/.accnRule

a0PN[i_] = a[i]/.accnRule/.omega0PNSqrule

(*velocity squared (v^2)/(c^2)*)
vSq=v[i] v[-i]/.velocityRule/.geometryRules

vSq=v[i] v[-i]/.velocityRule/.geometryRules/.omega0PNSqrule

vSqRule=(v^2)->vSq

Source  multipole  moments

(* Define Symmetric trace free projection of rank-2 tensor *)
DefTensor[STF2[i,j],Σ,Symmetric[{i,j}]];

(*Mass quadrupole (0PN)*)
Iij[-i,-j]:=m ν STF2[-i,-j] ( x[i] x[j] )

Print["Mass quadrupole ="];
Iij[-i,-j]

Derivative  of  multipole  moments

DefConstantSymbol[t];

(* Custom derivative rules *)
td[a_+b_]:=td[a]+td[b];

td[a_ b_]:=td[a] b+a td[b];

td[a_^n_]:=n a^(n-1) td[a];

td[G]:=0;
td[m]:=0;
td[ν]:=0;
td[γ]:=0;
td[r]:=0;
td[c]:=0; 
td[X1] :=0;
td[X2] :=0;
td[x_?NumericQ]:=0;

(* STF projector treated as constant *)
td[STF2[__]]:=0;

(* Orbital variables *)
td[x[i_]]:=v[i];
td[x[-i_]]:=v[-i];

td[v[i_]]:=a[i];
td[v[-i_]]:=a[-i];
td[a[i_]] := -(ω^2)v[i];

(* First derivative of Iij*)
Iijdot=td[Iij[-i,-j]]

(* Second derivative of Iij*)

Iijdotdot=td[Iijdot]

insideIijdotdot = Iijdotdot/(m ν STF2[-i,-j])

insideIijdotdot/.accnRule

insideIijdotdot/.accnRule/.omega0PNSqrule

insideIijdotdotN = insideIijdotdot/.accnRule/.omega0PNSqrule;

Iij2 = (m ν STF2[-i,-j])*insideIijdotdotN

(* Third derivative of Iij*)

Iijdotdotdot = td[Iij2]

insideIijdotdotdot = Iijdotdotdot/(m ν STF2[-i,-j])//Expand

insideIijdotdotdot/.accnRule

insideIijdotdotdot/.accnRule/.omega0PNSqrule /.stfSymRule1

(insideIijdotdotdot/.accnRule/.omega0PNSqrule //Simplify)/.stfSymRule1

(* STF Rule *)
stfSymRule1=v[j] x[i]+v[i] x[j]:>2 x[i] v[j]

Coefficient[insideIijdotdotdot/.accnRule/.omega0PNSqrule, -(((4 G m))/(r^3))] * -(((4 G m))/(r^3))

(Coefficient[insideIijdotdotdot/.accnRule/.omega0PNSqrule, -(((4 G m))/(r^3))] * -(((4 G m))/(r^3)))/.stfSymRule1

insideIijdotdotdotN = (Coefficient[insideIijdotdotdot/.accnRule/.omega0PNSqrule, -(((4 G m))/(r^3))] * -(((4 G m))/(r^3)))/.stfSymRule1;

Iij3 = (m ν STF2[-i,-j])*insideIijdotdotdotN

Energy  loss

Iij3*Iij3

STFcontractRule1=STF2[-i,-j] STF2[-i,-j]x[i] v[j]x[i] v[j]:>(1/2) (r^2)(v^2)

Iij3*Iij3/.STFcontractRule1

(Iij3*Iij3/.STFcontractRule1)/.vSqRule

Iij3Sq = (Iij3*Iij3/.STFcontractRule1)/.vSqRule;

dEdt0PN = ((-G)/(c^5))((1/5) Iij3Sq)